Optimizing Dynamic Antibiotic Treatment Strategies against Invasive Methicillin-Resistant Staphylococcus Aureus Infections using Causal Survival Forests and G-Formula on Statewide Electronic Health Record Data

Developing models for individualized, time-varying treatment optimization from observational data with large variable spaces, e.g., electronic health records (EHR), is problematic because of inherent, complex bias that can change over time. Traditional methods such as the g-formula are robust, but must identify critical subsets of variables due to combinatorial issues. Machine learning approaches such as causal survival forests have fewer constraints and can provide fine-tuned, individualized counterfactual predictions. In this study, we aimed to optimize time-varying antibiotic treatment –identifying treatment heterogeneity and conditional treatment effects– against invasive methicillin-resistant Staphylococcus Aureus (MRSA) infections, using statewide EHR data collected in Florida, USA. While many previous studies focused on measuring the effects of the first empiric treatment (i.e., usually vancomycin), our study focuses on dynamic sequential treatment changes, comparing possible vancomycin switches with other antibiotics at clinically relevant time points, e.g., after obtaining a bacterial culture and susceptibility testing. Our study population included adult individuals admitted to the hospital with invasive MRSA. We collected demographic, clinical, medication, and laboratory information from the EHR for these patients. Then, we followed three sequential antibiotic choices (i.e., their empiric treatment, subsequent directed treatment, and final sustaining treatment), evaluating 30-day mortality as the outcome. We applied both causal survival forests and g-formula using different clinical intervention policies. We found that switching from vancomycin to another antibiotic improved survival probability, yet there was a benefit from initiating vancomycin compared to not using it at any time point. These findings show consistency with the empiric choice of vancomycin before confirmation of MRSA and shed light on how to manage switches on course. In conclusion, this application of causal machine learning on EHR demonstrates utility in modeling dynamic, heterogeneous treatment effects that cannot be evaluated precisely using randomized clinical trials.


Introduction
Antimicrobial resistance (AMR) is a global public health threat.Every year, there are over 2.8 million antimicrobial-resistant infections in the United States, and at least 35,000 people die from these infections (Centers for Disease Control and Prevention, 2019).Identifying optimal antibiotic treatment regimens is key to increasing the chance of favorable clinical outcomes for individuals, and it can also have a role in preventing AMR at the population level.The World Health Organization (WHO) and the US Centers for Disease Control and Prevention (CDC) consider methicillin-resistant Staphylococcus aureus (MRSA) infection to be a bacterial disease of priority concern, because MRSA can affect multiple tissues and organs, and its treatment options are few, due to the fact that it is resistant to nearly all beta-lactam antibiotics (Centers for Disease Control and Prevention, 2019;World Health Organization, 2017;Lee et al., 2018;Turner et al., 2019;Rodvold and McConeghy, 2014;Alexander et al., 2023).
Invasive MRSA is a life-threatening complication in which the infection has spread inside the body, e.g., the bloodstream, deep-seated skin, or lungs (Lee et al., 2018).With invasive infections, patients often need emergent antibiotic treatment before the pathogen causing the infection is known.This initial treatment is called 'empiric treatment'.Clinicians use symptoms and source of infection to select empiric treatment according to experience and international guidelines (VanEperen and Segreti, 2016).Once the pathogen is identified by test results, clinicians can prescribe the first directed treatment.Currently, a bacterial identification culture test takes 2 -3 days.Because of this time delay from testing, in clinical practice, several sequential treatment scenarios are possible for each patient (Jeffres, 2017;Leekha et al., 2011).
Vancomycin is known as the first-line treatment for patients with a suspected MRSA infection (VanEperen and Segreti, 2016;Levine, 2006).Even though vancomycin has broad coverage, there can be multiple reasons to change vancomycin dose or switch to other antibiotics, including a lack of clinical improvement, updated information regarding antibiotic susceptibility, development of adverse effects, and even financial reasons (Jeffres, 2017).Nephrotoxicity is one of the main adverse events that patients experience while receiving vancomycin.When a patient experiences renal toxicity, doctors prescribe alternative antibiotics (Liu et al., 2020;Sawada et al., 2018;Filippone et al., 2017).
Although a randomized clinical trial (RCT) can be ideal approach for comparing first directed treatment for known drug susceptibility, this approach is cumbersome for dynamic, i.e., time-varying, antibiotic treatments especially in the context of delayed pathogen culture and antibiotic susceptibility testing.It is even more challenging to conduct a RCT with MRSA.MRSA must be treated immediately under multifaceted conditions.Estimating dynamic, heterogenous treatment effects among different patient groups is difficult to establish using a RCT, requiring a complex design and larger sample.However, observational data could be leveraged to calculate dynamic individualized treatment effects since numerous treatment scenarios are possible and observed in practice.
Previous studies of MRSA patients have compared vancomycin to other antibiotics, including daptomycin and linezolid (Schweizer et al., 2021;Yue et al., 2016;McCreary et al., 2020;Yeager et al., 2021;Paul et al., 2015;Liang et al., 2014).In a retrospective study of veterans with MRSA bloodstream infections, patients who were switched from vancomycin to daptomycin during the first three days after starting treatment had a lower risk of 30-day mortality than patients who did not switch from vancomycin (Schweizer et al., 2021).When compared to vancomycin, linezolid has shown promise for treating patients with skin and soft tissue infections (Yue et al., 2016) and ventilator-associated pneumonia (Peyrani and Ramirez, 2015).
Overall, there is a lack of studies that looked specifically at dynamic sequential treatment effects.To our knowledge, none applied methods for heterogeneous effect estimation.In this work, we analyzed longitudinal electronic health record (EHR) data from a large, statewide hospital settings, together with causal survival forest and g-formula, to estimate both average effects and heterogeneity of conditional effects for dynamic antibiotic therapy -vancomycin vs. others-on mortality in patients admitted with invasive MRSA infection.

Methods
First, we describe our ethical approvals, the data source, the study population, the longitudinal design, and the sequential treatment strategies, outcome, and covariates.Second, we give an overview of the causal survival forest and the g-formula methods, which we used to estimate dynamic treatment effects.

Ethics Statement, Data Source, and Derivation of Study Population
As authors, we abide to the ethical principles for medical research involving human subjects outlined by the World Medical Association in the Declaration of Helsinki.This study was reviewed and approved by the University of Florida's (UF) institutional review board (IRB) (protocol number IRB 201900652).We used deidentified data from a large university hospital system in Florida, UF Health, that comprises two primary hospitals in Gainesville and Jacksonville, as well as forty-five outpatient clinics in the state.Since 2011, UF Health uses the Epic system (https://www.epic.com/),and the EHR data is warehoused in the Integrated Data Repository (IDR, https://idr.ufhealth.org/).The IDR includes patients' demographics, clinical diagnoses, procedures, laboratory tests, and medications.Clinical diagnoses and procedures are encoded using the International Classification of Disease (ICD, https://www.who.int/standards/classifications/classification-of-diseases)ontology, 9th and 10the revision, while laboratory tests and medications are encoded via the Logical Observation Identifiers Names and Codes (LOINC, https://loinc.org/) and the RxNorm (https://www.nlm.nih.gov/research/umls/rxnorm/index.html)terminology, respectively.Data requests made to the IDR staff (https://idr.ufhealth.org/research-services/)should be in compliance with institutional, state and Federal regulations.The authors of this work are willing to share the study protocol and data analysis code.
Our study population includes adults (18 years and older) admitted to the hospital and diagnosed with invasive MRSA (the first one recorded).An MRSA diagnosis was confirmed with a culture test based on a biological sample, including blood, fluid, bone, kidney, liver, heart, lung, pancreas, etc.Individuals who had at least one-year of medical records before the identification of MRSA were included in this study to account for relevant medical history.Patients were followed during three sequential time points: (1) empiric treatment, (2) possible switch to the first directed treatment, and (3) the sustaining therapy with other switching options.The flowchart of the inclusion criteria for the study population is given in Figure 1.

Three Timepoints for Constructing Sequential Strategies
Selecting appropriate timepoints is crucial not only for assessing the sequential treatment effects, but also for identifying when treatment changes can be acted upon.We focused on the sequence of (1) empiric, (2) empiric to directed, and (3) sustaining treatment assignments as illustrated in Figure 2. "Time 1" refers to the interval period from admission to the receipt of culture test results.For example, if a patient received vancomycin during this interval, we labeled the patients into vancomycin group at Time 1.During this time, the definitive organism and antibiotic susceptibility test results are not known.We collected the relevant measurement proxies from the EHR to ascertain the treatment propensity, as well as potential causes of early/late adverse reactions that can entail contraindications."Time 2" is a measure of the preliminary response to the empiric antibiotic therapy, also known as early response assessment.Since the initial response period is typically assessed within 3 -7 days, we fixed it at 3 days after the culture test.With these results, providers may continue with their empiric treatment prescription (perhaps with dose adjustment) or switch to the directed treatment.Various clinical factors, such as nephrotoxicity, may affect this transition."Time 3" involves monitoring the antibiotic treatment and sustaining therapy for the recommended time.For MRSA, this is typically between 7 and 14 days after the initiation of therapy (or even longer) depending on the severity and location of the infection.During this time period, the overall effectiveness of the therapy, any remaining signs of infection, and the potential for recurrence or complications is assessed.In our study, we define this third time point as the 7 days from the first directed treatment (i.e., "Time 2" + 7 days, or "Time 1" + 10 days).
At each of the three timepoints, we assessed if the patient was taking vancomycin or they were prescribed another.In total, 8 different sequential treatment strategies were considered.For example, if a patient started with vancomycin as empiric treatment and maintained the same treatment at timepoint 2, but then changed to another antibiotic at timepoint 3, the patient would have '1-1-0' as the value for the sequential treatment strategy variable.

Study Intervention, Outcome, and Covariates
We defined two different interventions that could be applied in clinical practice: one was a three-point treatment sequence (modelled using a causal survival forest) and the other was a treatment update at a given time point (modelled through g-formula).For the first intervention, we encoded a binary treatment variable to indicate whether there was an antibiotic change between the previous time point to next time point (i.e., treatment change from Time 1 to Time 2 or treatment change from Time 2 to Time 3, and any change between the first and another time point, i.e., from Time 1 to Time 2/Time 3).For the second intervention, vancomycin was the target treatment, and any other antibiotic was pooled into the control group, corresponding to the 8 sequential treatments.The study outcome was the time from the onset of bacterial infection (i.e., culture collection date set as the index date) to death or discharge within a 30-day horizon (i.e., 30-day mortality).Study covariates were both time-fixed for the first intervention and time-varying for the second intervention.Time-fixed covariates measured before MRSA onset or at index date included patient's demographics (age, sex, race), Charlson's comorbidity index, admission type, intensive care unit (ICU) stay, healthcare acquired infection, and previous antimicrobial resistance testing.We also collated all prior clinical diagnoses present with at least 10% frequency in the study population (mapping all ICD-10 codes into ICD-9), to investigate additional potential drivers of the outcome, as done in another study (Jun et al., 2022).The setup for time-varying covariates collected after the index date (time-fixed) is illustrated in Figure 3 using a causal directed acyclic graph (DAG).For example, one time-varying confounder is the nephrotoxicity variable, defined as a 50% decrease in creatinine clearance (CrCl) from a baseline value (Wong-Beringer et al., 2011).If a patient was missing the creatinine clearance value at a given time point, the previous creatinine value closest to the time point was used.

Causal Survival Forests and G-formula
Causal survival forests (CSF) are an adaptation of the causal forest algorithm, a nonparametric method for estimating heterogeneous treatment effects in survival settings with right-censored data (Cui et al., 2023).The causal effect of the antibiotic sequential strategy is estimated under the Neyman-Rubin's potential outcome model framework (Imbens and Rubin, 2015).This estimation operates within a statistical setting where we have n independent and identically distributed subjects (i = 1, … , n), and we observe each subject's tuple (X i , Y i , W i , D i ).
In this context, X i denotes a vector of covariates, Y i represents the observed response variable (i.e., days to death), W i ∈ {0, 1} signifies the binary treatment assignment (specifically, whether or not a change in the antibiotic was made), and D i acts as an event indicator, signifying whether the event (i.e., death) took place.Given this configuration, we can identify (1) if there are heterogeneous effects among patients and (2) which specific population (i.e., which combination of covariates) shows high heterogeneity.We have employed the concepts of the Conditional Average Treatment Effect (CATE) and the Rank-Weighted Average Treatment Effect (RATE) to quantify these heterogeneous effects.
The CATE is defined by the equation and is the expected mean of difference between potential outcomes Y i (0), Y i (1) given auxiliary covariates X i .
CATE can be used to derive treatment prioritization rules, and the RATE serves the purpose of evaluating how good treatment prioritization rules are at distinguishing sub-populations with different treatment effects, or whether there exists notable heterogeneity.RATE only considers ranking of each patient's rather than considering numeric size of the score.For quantifying the treatment benefit, the Targeting Operator Characteristics (TOC) is further calculated.
In the TOC equation, TOC(q) implies the top q-th fraction of individuals with the largest prioritzation score S(X i ).F S(Xi) is the distribution function of S(X i ) for comparing the ATE in the top q-th fraction of individuals with the largest prioritization score S(X i ) with the overall ATE from treating everyone.If the TOC is equal to 0, it means that there is no benefit in stratifying the treatments using given prioritization rules.The parametric gformula is an extended version of standardization by Robins for time-varying treatments and confounders.It uses the identification assumptions of inverse probability weighting, but it models the outcome means instead of the treatment equation (Ezzati et al., 2004;Westreich et al., 2012;McGrath et al., 2020).All analyses were conducted using the R software (https://www.r-project.org/),including the 'grf' (Tibshirani et al., 2023) and 'gfoRmula' (McGrath et al., 2020) packages.

Population Characteristics, Outcomes, and Dynamic Treatment Assignments
Among 1,433 patients admitted between 2011 and 2019, with a confirmed MRSA diagnosis, 914 had at least one year of prior medical history from the onset of the infection and 872 patients had complete record data from admission to discharge.Of these, Time 1 was observed in 817 patients, Time 2 in 707 patients, and Time 3 in 427 patients.As the objective was to reach the sustained treatment time point, the final study population comprises the last subset of 427 patients.
The mean age of the study population was 55 years, 48.9% were male, 60.4% white, and 34.4% had multi-drug resistance (more than three antibiotic classes).The proportion of patients who stayed in the ICU was 52.2% and 24.1% were assumed to have an healthcareacquired infection.The overall length of admission was median 19 days, from a minimum of 10 days to a maximum of 1,050 days.Out of the 427 patients, 33 patients (7.7%) died within 30 days from the MRSA onset.The summary statistics on the study population are given in Table 1.
Among the 427 patients, 96.9% started with vancomycin, 56.2% used vancomycin treatment throughout all three timepoints, while 43.8% of them changed from vancomycin to another antibiotic at least once; 21.3% changed their treatment from Time 1 to Time 2 and 32.6% from Time2 to Time 3. In Table 2 we summarized the number of subjects within each possible scenarios of antibiotic treatment at the three sequential timepoints corresponding to the empirical, directed, and sustaining treatment periods, along with the proportions of suspected nephrotoxicity, which is one of the main reasons for antibiotic change.

Estimation of Treatment Effects
For the first intervention, we generated two estimates of the conditional average treatment effect using the CSF.Model 1 included the expert-selected variables of the DAG (sex, race, age, Charlson's comorbidity score, ICU stay, and multidrug resistance).The average treatment effect for changing antibiotics during any timepoint showed a reduction in the probability of death (Mean = −0.07619,SE=0.02953).All covariates had no effect on mortality (i.e., p-value below 0.05).Model 2 also included the DAG variables, but expanded all individual comorbidities.The average treatment effect of this model showed a negative effect (Mean = −0.0810,SE=0.02953) similar to Model 1. Table 3 details CATE results using the best linear predictor for both Model 1 and Model 2.
We then assessed heterogeneity for the CSF models using the area under the TOC curve (AUTOC) from the RATE function (Figure 4).In Model 1, there was minimal heterogeneity of treatment effect (AUTOC = 0.01, S.E = 0.06).In Model 2, the group with quantile above 15 would benefit from the treatment change (AUTOC = −0.04,S.E = 0.06).
In the analysis of the second intervention, we used the g-formula to estimate the effect of vancomycin vs. other treatment updated at each time point on 30-days mortality, in the presence of time-varying and time-fixed confounding.Table 4 shows the g-formula mean, ratio, and difference for the reference Intervention 0, which was the observed treatment (natural course), compared to the Intervention 1 scenario of "Never treat with vancomycin or switch immediately after empiric treatment", and Intervention 2 scenario of "Always treat with vancomycin and do not switch to another antibiotic".The mean ratio and mean difference highlight the difference between intervention scenarios using observed scenario as a reference group.For the whole study population, under the reference scenario, the empirical risk of the outcome was 7.7%, corresponding to the event proportion (30-day mortality) in our study population.If all patients were assigned to Intervention 1, the risk increased approximately to 14.6%.If all patients were assigned to Intervention 2, the risk was lower than the observed treatment scenario which was about 5.0%.
In the subset of patients who started their empiric treatment with vancomycin, under Intervention 1, the risk increased to 16.7% from the reference scenario.Under Intervention 2, the risk was 4.97%, lower than the observed treatment scenario.Overall, results in the whole and vancomycin-empiric population were similar, and the current sequential vancomycin treatment was better than not giving anyone vancomycin, but worse than giving everyone vancomycin.

Discussion
We found that switching from vancomycin to another antibiotic improved survival probability.Additionally, there was benefit from initiating vancomycin compared to not using it at any time point.Our findings are consistent with the general knowledge of clinical efficacy of vancomycin in the treatment of invasive MRSA infections obtained through RCTs.In our population, over 95% of patients were prescribed empiric vancomycin treatment, reflecting the common preference for this antibiotic in managing potential infections that have not yet been confirmed (VanEperen and Segreti, 2016).However, despite being the first-line choice of treating invasive MRSA infections, vancomycin was consistently used during all three timepoints in only half (56%) of the of the population and many patients who initially treated with vancomycin switched to another antibiotic at least once.This change might suggest concerns over vancomycin's side effects, particularly nephrotoxicity, from the provider (Jeffres, 2017), reflecting the importance of considering these time-varying components in our models.This is also supported by the increased prevalence of nephrotoxicity in those who received consistent vancomycin treatment.Despite this, we found a decreased mortality probability when patients stayed on vancomycin, emphasizing its role despite potential complications.
In addition to time-varying treatment effects, we investigated putative treatment heterogeneity, which is key to the development of personalized treatments that cater to individual patient characteristics (Varadhan et al., 2013).We found low evidence of heterogeneity, although it could have been due to lack of power.Lack of power is even more problematic with RCTs, that can become cumbersome and resource-heavy to include sufficiently diverse populations.Furthermore, RCTs can include inherent selection bias because of strict inclusion criteria.For example, those with establish renal failure are typically excluded from vancomycin trials (Paul et al., 2015), however these individuals remain a priority group in assessing effects of vancomycin due to known side effects of the antibiotic.
Our study has a number of limitations.Firstly, our study design and data analysis make simplifications (albeit clinically reasonable) in the treatment staging and decision-making process that are a conceptual abstraction.The timepoints utilized in this study do not exactly align with real clinical settings and patient populations.While culturing methods remain the preferred method for confirming MRSA infections, different institutions may implement other approaches that would deviate from the timeline we defined.It is important to note here that, even if the measurement time points are the same across all patients, there can be also chance of including immortal time bias which differs among treatment paths.Secondly, our analysis does not consider the dosage of vancomycin therapy.This is particularly important to take into account in future analyses, due to its relationship with nephrotoxicity that influence on antibiotic selection.Thirdly, we considered all MRSA infections in this analysis and did not differentiate between specific types of infections (e.g., bloodstream, lung, skin and soft tissue).It is possible that the relationship between sequential vancomycin therapy and outcomes may vary between anatomic sources of infections due to the inherent pharmacodynamics of vancomycin.
As statistical methodologies continue to advance, promising opportunities arise for future research in sequential treatment optimization, particularly through the incorporation of algorithms from causal inference and machine learning.Examples include the flexible Bayesian Additive Regression Trees (BART) and Counterfactual Regression (CFR), which adaptly handle high-dimensional environments and intricate non-linear relationships, respectively (Hill, 2011;Shalit et al., 2017).Despite not being originally designed for time-varying treatments or confounders, innovative variations of these models have been developed (Linero and Zhang, 2022).In our study, we focused on sequential time-varying treatment options and confounders, utilizing the G-formula due to its ability to manage complex, dynamic scenarios.However, evaluating a broader array of models could further substantiate our research findings.

Conclusion
In this study we operationalized sequential treatment strategies aimed at identifying relevant heterogeneity and optimizing risk based on individual patients' characteristics.We demonstrated the utility of applying causal machine learning to real-world data within a framework that can be used to screen multiple intervention hypotheses, especially for life-threatening conditions, and select the most promising to be tested with conventional RCTs, possibly saving resources and lives.Description of treatment timeline and three timepoints Note: T1 is an interval period from the start of admission to the culture test result confirmed, T2 is the 1st day of the directed treatment, and T3 is the 8th day of directed treatment    Number of study subjects within each possible scenarios of antibiotic treatment (A 1 A 2 A 3 ) at the three sequential timepoints corresponding to the empirical, directed, and sustaining treatment periods (T 1 T 2 T 3 ), together with the proportion of suspected nephrotoxicity that can trigger antibiotic change.

Figure 1 :
Figure 1: Flowchart of inclusion criteria to derive the study population Figure 2:

Figure 3 :
Figure 3: The directed acyclic graph representing causal relationships among sequential antibiotic treatment (A1 at Time 1, A2 at Time 2, and A3 at Time 3), influencers and confounding factors, and 30-day mortality (Y) Note: A1, A2, A3 are treatment variables at each time point (either vancomycin or others).Y is a 30-day mortality.Demographics, comorbidity, severity, and past multi-drug resistance from previous EHR records are measured at baseline.Nephrotoxicity is a time-varying confounder that affects the following antibiotic prescriptions (A2, A3).In this figure, we included known risk factors between antibiotic treatment and mortality.

Figure 4 :
Figure 4: Targeting Operator Characteristics Curves for CSF-based CATE models: Model1 depicted on the left side and Model2 depicted on the right side graph

Table 1 :
Baseline characteristics of the vancomycin started groups and overall population

Table 4 :
G-formula estimation the effect of vancomycin vs. other treatment updated at each time point on 30-days mortality, in the presence of time-varying and time-fixed confounding, stratified by empiric treatment (Time 1).